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ABSTRACT 

Relativistic collisionless shocks in electron-ion plasma are thought to occur in the afterglow phase 
of Gamma-Ray Bursts (GRBs), and in other environments where relativistic flows interact with the 
interstellar medium. A particular regime of shocks in an unmagnetized plasma has generated much 
interest for GRB applications. In this paper we present ab-initio particle-in-cell simulations of unmag- 
netized relativistic electron-ion shocks. Using long-term 2.5-dimensional simulations with ion-electron 
mass ratios from 16 to 1000 we resolve the shock formation and reach a steady-state shock structure 
beyond the initial transient. We find that even at high ion-electron mass ratios initially unmagnetized 
shocks can be effectively mediated by the ion Weibel instability with a typical shock thickness of ~ 50 
ion skin-depths. Upstream of the shock the interaction with merging ion current filaments heats the 
electron component, so that the post shock flow achieves near equipartition between the ions and elec- 
trons, with electron temperature reaching 50% of the ion temperature. This energy exchange helps 
to explain the large electron energy fraction inferred from GRB afterglow observations. 

Subject headings: acceleration of particles - collisionless shocks - gamma rays: bursts 



1. INTRODUCTION 

Relativistic collisionless shocks propagating in unmag- 
netized electron-ion plasmas are an essential ingredient 
in the theory of GRB afterglows (see Piran 2005 and 
Waxman 2006 for reviews). These shocks are expected 
to generate sub-equipartition magnetic fields and accel- 
erate nonthermal particles that are responsible for the 
observed synchrotron emission. Weibel instability has 
emerged as the leading mechanism for shock formation in 
weakly magnetized plasmas (Medvedev and Loeb 1999; 
Gruzinov and Waxman 1999). This instability can gener- 
ate small scale magnetic fields in counterstreaming plas- 
mas and can isotropize the flow. Weibel instability was 
simulated by a number of authors in various contexts 
(e.g., Nishikawaet al. 2003; Silvaet al. 2003; Frederiksen 
et al. 2004; Medvedev et al. 2005; Nishikawaet al. 2005). 
Although all groups observe the initial filamentation, suf- 
ficiently large simulations that lead to shock formation 
have only been done in the case of electron-positron pair 
shocks (Spitkovsky 2005, further S05; Chang et al. 2007, 
further CSA07). For the electron- ion plasmas the pio- 
neering simulations of Frederiksen et al. (2004), Hededal 
et al. (2004) and Hededal (2005) were not large enough 
to see the complete thermalization of the ions and the 
shock formation. This led Lyubarsky & Eichler (2006) to 
question whether the Weibel instability in the electron- 
ion plasma is fast enough to generate a shock transition 
that is thinner than the Larmor radius in the weak ISM 
field. In this paper we demonstrate via large numeri- 
cal simulations that the ion Weibel instability is indeed 
very effective at establishing the shock transition in an 
unmagnetized electron-ion plasma even for realistic mass 
ratios. We study the shock structure and energetics of 
the flow and find that electron heating in the foreshock 
is an important catalyst to shock formation via filamen- 
tation instability. From simulations we find the energy 
fraction that is transferred to the post-shock electrons 
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and compare it to the values inferred from GRB afterglow 
observations. In ^2] we describe the particle-in-cell sim- 
ulations of shock formation, and then discuss the shock 
structure in ^3] and electron heating in Q 

2. SIMULATION SETUP 

Particle-in-cell (PIC) simulations have proven to be 
a valuable tool in understanding collisionless plasma 
physics from first principles (Birdsall and Langdon 1992). 
For our study we use a relativistic electromagnetic PIC 
code TRISTAN-MP, which is a parallel descendant of 
the code TRISTAN (Buneman 1993). Among the mod- 
ifications we introduced are improved noise properties 
of the charge deposition and the suppression of numeri- 
cal Cherenkov instabilities of ultrarelativistic particles. 
This enables the propagation of cold relativistic plas- 
mas through the grid for thousands of plasma times 
(Spitkovsky and Arons 2007, in prep, hereafter SA07). 
To simulate shock formation we reflect a relativistically 
moving cold electron-ion stream from a conducting wall. 
This is equivalent to colliding two streams of identical 
plasma head-on, but saves 1/2 of the computational ef- 
fort. We validated this setup by comparing with a true 
head-on collision. In hydrodynamics a collision of two 
plasmas should create a contact discontinuity and two 
shock waves propagating away from the contact into the 
"upstream." The wall in our simulation represents the 
contact, and we simulate only a moving forward shock. 
Our simulation is performed in the "downstream" frame 
of the contact, not in the shock frame. Since the down- 
stream plasma is at rest, no extra boosts are needed 
to study the spectra. This helps to disentangle a rel- 
ativistically moving system of precursors, shocks and 
contacts that forms in collisions of jets with station- 
ary plasmas (Frederiksen et al. 2004, Nishikawa et al. 
2005). The characteristic scales of interest in an electron- 
ion shock are on the order of tens of ion skin depths 
(c/ujpi = (47re 2 n/7miC 2 ) _1//2 , for a plasma with density 
n and characteristic ion energy jrriiC 2 ). We typically 
use 10 cells per electron skin depth, so in order to re- 
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Fig. 1. — Density profiles during the formation of rrii/m e = 500 
shock. The incoming flow is moving to the left with 70 = 15. 

solve the ion scales for realistic mass ratios we would 
require simulations with several thousand cells on the 
side in 3D and tens of billions of particles. This re- 
quires a dedicated large-scale simulation effort. Instead, 
we have chosen to simulate shock formation in 2D, en- 
couraged by the similarity between 2D and 3D shocks 
in pair plasma (CSA07, S05). The simulation plane in- 
cludes the direction of streaming and one transverse di- 
mension. The magnetic field is out of the plane (CSA07). 
We simulated shocks in an initially cold, unmagnetized 
plasma moving with Lorentz factor 70 = 15 for a range of 
ion-electron mass ratios: rrii/m e = 16, 30, 100, 500, 1000. 
The transverse box size varied so that the box encom- 
passed ~ 40 upstream ion skin depths (2048 cells for 
rrii/me = 16,30; 4096 cells for mi/m e = 100; 8192 cells 
for rrii/m e = 500, 1000). The longitudinal size expanded 
with the simulation time up to 10 5 cells for rrii/m e = 100, 
corresponding to 10 3 c/uj p i (this was our longest run). We 
typically used 2 particles per cell per species in the unper- 
turbed upstream, which increased to 6 behind the shock 
due to shock compression. We have tested a larger num- 
ber of particles per cell and larger transverse box sizes, 
but the results did not appreciably change. 

3. SHOCK STRUCTURE AND EVOLUTION 

In fig. [I] we display the time sequence of transversely 
averaged density profiles that show the formation of a 
shock in a rrii/m e — 500 plasma (the same occurs for 
other mass ratios). Initially the flow with density n\ 
and 70 = 15 is moving to the left. After the bounce 
from the wall the reflected and the incoming plasmas 
stream through each other, increasing the density to 
2m (fig. [lji) . The electrons undergo Weibel insta- 
bility and tnermalize to their upstream kinetic energy, 
but the ions are still cold. After ~ lOOccC 1 the den- 
sity near the wall begins to rise, as the ions are ran- 
domizing (fig. [l] b). The ions at the head of the re- 
flected flow are still cold, and propagate close to c into 
the upstream. This is the "initialization precursor" - 
fast particles that are the remnants of the initial col- 
lision. They appear as a moving density bump that 
always outruns the shock (fig. [lb-f ). With time this 
bump is eroded as the particles aeceler ate. We define 
the shock as the density compression that propagates 
away from the wall in fig. [l]>f. The shock satisfies the 
hydrodynamic jump conditions after the initial transient: 
n 2 /n 1 = T ad /(T ad - 1) + l/(7o(r ad - 1)) = 3.13, where 
the adiabatic index is Y ad = 3/2, appropriate for a 2D 
relativistic gas. 

When the initialization precursor erodes so that the 
density in front of the shock is close to the unperturbed 
upstream density, the shock becomes steady. In this 
stage the integrated quantities through the shock do not 
significantly evolve as the shock moves though the grid. 
The shock structure snapshot in this stage is shown in 




Fig. 2. — Steady state structure of mi/m e = 100 collisionless 
shock, a) Density structure in simulation plane (normalized by the 
unperturbed upstream density); b) Magnetic energy density cb, 
in units of the upstream energy density; c) Transversely averaged 
plasma density; d) Transversely averaged cb- 

fig. [2] for a 250c/ (j p i section around the shock from the 
rrii/m e = 100 simulation. The shock speed is close to 

Vshock = c(T ad - 1)^(70 - l)7(7o + 1) = 0.47c, in accord 
with jump conditions. Plasma density (fig. [2^1) shows fil- 
amentation in the upstream region with filaments reach- 
ing lOc/cjpi in size. Magnetic energy (J2Jd) is also filla- 
mentary with enhancements at the edges of density fila- 
ments. Filamentation is driven mainly by the ion dynam- 
ics, and the shock compression starts where the filaments 
merge and disintegrate. The shock transition, which is 
~ 50c/uj p i thick (fig. 2 2) corresponds to a peak in the 
magnetic energy (fig. 2i). Incoming flow is isotropized 
by this magnetic field. vVhile the magnetic energy den- 
sity in the shock can reach local equipartition with the 
upstream flow energy (e# = £> 2 /47T7oni m^c 2 ~ 1), the 
average magnetic energy is near 10% — 15% of equipar- 
tition at the shock. Magnetic filaments are destroyed 
in the shock, and the downstream magnetic field forms 
islands as in the pair shock simulations (CSA07). The 
field energy decays below 4 x 10 -3 of the upstream en- 
ergy density; however, runs with more particles/cell are 
needed to reliably study the downstream field further. 

The filaments of density and magnetic field are not 
stationary in front of the shock. They are advected with 
the upstream flow and merge on the transit time of a 
filament towards the shock. As a result, the individual 
clumps that enter the shock change over time; however, 
the average density profile is relatively stable. 

4. ELECTRON HEATING 

In the steady state the shock is not influenced by the 
wall or the initial collision. It is a self-propagating struc- 
ture. In order to maintain continuing filamentation in 
the upstream the cold incoming fluid should experience 
counterstreaming, which is provided by a population of 
particles escaping from the downstream (Milosavljevic 
et al. 2006, SA07). In fig. [3] we show the longitudinal 
momentum space for ions (fig. ^jp) and electrons (fig. 
[3J3) for rrii/m e = 100 run. The incoming flow has nega- 
tive values of four- velocity. In this snapshot the flow is 
stopped and thermalized at the shock for x < 250c/uj p i 
(fig. |3k shows the density structure of the shock for ref- 
erence). There is a clear population of particles with 
positive four-velocities streaming away from the shock. 
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These particles, though being hotter and more tenuous 
than the incoming fluid, cause the filamentation instabil- 
ity in the upstream. A structure with positive 4- velocity 
at x > 450c/ 'cOpi in fig. [3]o is the remnant of the ini- 
tialization precursor. Over time its contribution to the 
initiation of filamentation diminishes (we have seen it 
disappear completely in runs with smaller mass ratios). 
The particles reflected from the shock occupy a region 
that extends > 300c/ 'co P i before the shock for the times 
to which we have evolved the simulation. 

From fig. [3]3 it is clear that the electrons thermalize 
with 7- factors that significantly exceed 70, suggesting 
electron heating before the shock. To quantify this ef- 
fect we plot the average energy per particle normalized 
by the upstream ion energy in fig. with red and blue 
curves for ions and electrons respectively. The solid lines 
show only the electrons and ions that are moving towards 
the shock; these lines show the partition of energy in a 
fluid element as it approaches the shock. The electrons 
gain 35% of the initial ion energy by the time they reach 
the downstream. In the units of downstream ion energy, 
the electrons are at 50% of equipartition with the ions. 
If we include the reflected particles into the mean en- 
ergy calculation (dashed line in fig. |3]i), both species 
show reheating near the shock, indicating that the re- 
flected particles are a separate population that get their 
energy by bouncing off the shock front and can spend 
more time near the shock than the transit time of the 
incoming fluid. The energy spectra of slices from the up- 
stream (fig. [3]p) and the downstream (fig. |3^) show that 
while in the upstream the ions were a cold beam and 
electrons were thermalized to their initial energy, in the 
downstream both species are near Maxwellian distribu- 
tions with similar temperatures. A non-Maxwellian tail 
can also be observed, especially in electrons, and is due 
to reflected particles (Spitkovsky 07, in prep). 

We studied the mechanism of electron heating by plot- 
ting the orbits of particles from simulations. The ion 
filaments are not completely neutral, with an excess of 
positive charge in the middle of the filament and a shield- 
ing negative charge on the periphery. Hence, there is an 
electric field which accelerates electrons to the center of 
the filament (in addition to giving the ExB motion with 
the filament towards the shock). In general, this accel- 
eration is reversible as the electrons lose energy on the 
way out of the filament (Hededal et al. 2004, Medvedev 
2006). However, the development of Weibel instability 
increases the charge of the filaments with time, so on 
every return particles experience a deeper potential well 
and increase the amplitude of energy oscillations. In the 
nonlinear stage of the instability (near the shock), the 
filaments merge and break up on the time scale of elec- 
tron oscillation in the filament. In this region the hot 
electrons can switch ion filaments while retaining their 
energy. Electrons are also attracted to isolated clumps 
of positive charge that result from the break up of ion 
filaments near the shock. When averaged in the trans- 
verse direction these clumps of charge form an attractive 
electrostatic potential well extending ~ 200c/ co P i in front 
of the shock that accelerates electrons and decelerates 
ions. The potential is variable on the time scale of tran- 
sit of charge clumps, so the electrons can retain some of 
the gained energy after exiting the well and crossing the 
shock. A detailed theory of the heating process will be 
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Fig. 3. — Internal structure of mj/m e = 100 shock, a) Average 
density profile; b) Momentum space density N(x,jf3 x ) for ions; c) 
Momentum space density for electrons; d) Average particle energy 
for particles that move towards the shock (in units of upstream ion 
energy): red - ions, blue - electrons. Dashed lines show average 
energy that includes reflected particles; e) Particle spectrum in the 
downstream slice centered on x = 150c/ ui p i, red for ions, blue for 
electons; f) Spectrum in the upstream slice at x = 530c/ 'uj p i. 

presented elsewhere (Medvedev & Spitkovsky, in prep). 

5. DISCUSSION 

In this paper we studied the structure of collisionless 
shocks in an initially unmagnetized relativistic electron- 
ion plasma. Our 2D PIC simulations are sufficiently 
large to resolve the formation of a shock as a density 
jump of thickness ~ 50c/ w P i which satisfies the hydrody- 
namic jump conditions. Self-generated transverse mag- 
netic fields reach 15% of the initial ion energy in the shock 
transition. An important ingredient is the foreshock re- 
gion that extends several hundred c/uj vi into the up- 
stream. Here, the incoming cold flow interacts with the 
hot tenuous plasma that escapes from the downstream. 
This interaction leads to the ion Weibel instability and 
the formation of filaments of current and density. Elec- 
trons oscillating in the electromagnetic field of the grow- 
ing and merging ion filaments are efficiently heated close 
to equipartition with the ions. This heating is essential to 
the formation of unmagnetized shocks. As was pointed 
out by Lyubarsky & Eichler (2006), ion Weibel instability 
slows down if electrons have temperature much smaller 
than the ion energy. Such electrons can efficiently screen 
the ion filaments from mutual attraction. Electron heat- 
ing reduces this shielding, increasing the electron Larmor 
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radius and skin depth, and allowing them to leave the fil- 
aments more easily. This increases the charge separation 
inside the filaments and facilitates merging. Electron 
acceleration in ion filaments was previously discussed by 
Hededal et al. (2004) and Medvedev (2006). In the shock 
context we can also point out the effects that can allow 
electrons to keep the gained energy: the evolution and 
merging of filaments on the transit time to the shock 
and the escape of energized electrons from the filaments. 
This creates a time- varying potential well in front of the 
shock, which contributes to net heating. 

The amount of heating we find is considerable - the 
downstream electron energy fraction is e e ~ 0.5. This 
supports the canonical picture of GRB afterglows, where 
electron energy fraction in the radiating electrons e e ~ 
0.1 is commonly inferred (Piran 2005). We note that the 
bulk of our downstream electrons is thermal (fig. [3]f), 
so it is reasonable to assume that the energy fraction in 
the nonthermal component, once it appears, would be 
smaller than what we obtain. Our result also constrains 
the minimum electron Lorentz factor in an accelerated 
power-law to be > jom p /m e (cf. Eichler & Waxman 
2005). The heating mechanism is robust even in the ab- 
sence of a shock: we performed 2D periodic box sim- 
ulations of counterstreaming electron-ion plasmas with 
7o = 15 for mass ratios of 100 and 500. In both cases 
the electrons gained 24% of the initial flow energy, which 
is smaller than in the shock setup, but still much larger 
than the original electron energy. 

We studied the properties of the shocks as a function of 
electron-ion mass ratio. Interestingly, for rrii/m e > 30, 
the properties such as shock width (in c/uo pi ) and mag- 
netic and electron energy fractions do not significantly 
change with ion mass. Although the largest mass ra- 
tio we tried was 1000, the convergence of shock proper- 
ties with mass suggests to expect no surprises even at 
rrii/m e = 1836. Also, shocks with 10 < 70 < 100 show 
very similar behavior, so our results for 70 = 15 are rep- 
resentative of ultra-relativistic shocks. 

The shock structure seems to have reached a steady 
state in our simulations. This means that we do not 
see a significant evolution in the integrated quantities 
as the shock propagates through the grid on time scales 
of our runs (up to lO 3 ^;" 1 ). However, there are some 
features that continue to evolve. There are particles that 
are escaping into the upstream at ~ c, and with time 
the spatial extent of the "cloud" of reflected particles 
increases. If there are Fermi-type processes in the shock, 
they would be associated with these particles (Spitkovsky 
07, in prep). Based on our simulations we cannot rule 
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out the formation of high energy particle population that 
could change the shock structure on longer timescales. 

A lingering question is whether our 2D simulations cap- 
ture the relevant 3D physics. To check this we performed 
a rrii/m e = 16 simulation of a 3D relativistic shock on 
a relatively small 256 2 x 5000 grid (only 6 2 x 12bc/u pi ). 
We obtain the density ramp and the formation of a shock 
with electron heating at 30% of the upstream energy that 
is very similar to the 2D case. From this we expect that 
our 2D results for higher mass ratios will hold in 3D. 

Electron heating to near-equipartition with the ions 
in the shock implies that relativistic electron-ion shocks 
should be similar to the electron-positron pair shocks, 
especially in the downstream region. There, both elec- 
trons and ions have similar effective mass and Larmor 
radii. The behavior of the downstream magnetic turbu- 
lence should then be very similar between the pair and 
electron- ion shocks. The magnetic islands in the down- 
stream of fig. |2Jd are on the scale 10c/ 'u P i (which is also 
close to 10c/a;p e in the downstream electron skin depths). 
Similar islands were observed in CSA07 for the down- 
stream of pair shocks. We expect that higher resolution 
simulations of the downstream region in e-ion shocks will 
find similar decay of the field with time as in CSA07. 

Recently, the interest in understanding the electron-ion 
temperature ratio in collisionless shocks has increased 
due to the new observations of T e /Ti in supernova 
remnant shocks (Rakowski 2006) and in cluster shocks 
(Markevitch & Vikhlinin 2007). The SNR observations 
imply 1 > T e /Ti > 0.1, decreasing with shock velocity 
(Mach number). Although we find T e /Ti w 0.5 in the 
relativistic case and no significant variation in our re- 
sults with shock Lorentz factor, we caution against direct 
comparison between our results and the nonrelativistic 
shocks. First, the SNR shocks may be magnetized, and 
here we only deal with unmagnetized shocks. Second, 
the nonrelativistic problem is more complicated because 
of the separation of scales implied by the different ther- 
mal velocities of ions and electrons. Other instabilities 
(such as electrostatic two-stream, firehose, etc) may also 
be important in nonrelativistic shocks. So, it is not sur- 
prising that the structure of these shocks may depend on 
more parameters. The nonrelativistic regime can, how- 
ever, be explored with the existing PIC technology. 
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